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We theoretically investigate the current-voltage {I-V) property of two-dimensional Coulomb 
blockade (CB) arrays by conducting Monte Carlo simulations. The I-V property can be di- 
vided into three regions and we report the dependence of the aspect ratio 5 (namely, the 
lateral size Ny over the longitudinal one Nx). We show that the average CB threshold obeys a 
power-law decay as a function of 5. Its exponent 7 corresponds to a sensitivity of the thresh- 
old depending on 5, and is inversely proportional to A^^. (i.e., 5 at fixed Ny). Further, the 
power-law exponent C, characterizing the nonlinearity of the I-V property in the intermedi- 
ate region, logarithmically increases as 5 increases. Our simulations describe the experimen- 
tal result C = 2.25 obtained by Parthasarathy et al. [Phys. Rev. Lett. 87 (2001) 186807]. In 
addition, the asymptotic I-V property of one-dimensional arrays obtained by Bascones et 
al. [Phys. Rev. B. 77 (2008) 245422] is applied to two-dimensional arrays. The asymptotic 
equation converges to the Ohm's law at the large voltage limit, and the combined tunneling- 
resistance is inversely proportional to 5. The extended asymptotic equation with the first-order 
perturbation well describes the experimental result obtained by Kurdak et al. [Phys. Rev. B 57 
(1998) R6842]. Based on our asymptotic equation, we can estimate physical values that it is 
hard to obtain experimentally. 
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1. Introduction 

A Coulomb blockade (CB)^'^^ emerges in condensed matter physics, and it causes thresh- 
old and nonlinear current-voltage {I-V) behavior. In some sense, CB can be regarded as 
a phenomenon that occurs in disordered systems with thresholds such as charge-density 
waves^^ and Wigner crystals.'^^ CB was first studied in a single-electron transistor. Nowa- 
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days, it is studied in many systems such as arrays of metallic islands, ^ metal nanocrys- 
tal arrays,^"^^^ molecular arrays, ^^"^^^ a Tomonaga-Luttinger liquid such as a carbon nan- 
otube,^^^^^ and graphene quantum-dot arrays. Some numerical studies of CB have been 
done, e.g. Monte Carlo (MC) methods,^*'"^"*^ molecular dynamics (MD) simulations,^^"^^^ and 
circuit dynamics.^^'^^^ In most of these cases, a core work is to investigate its I-V property. It 
is known that the I-V property of typical CB arrays can be divided into three characteristic 
regions according to path flow of electrons: near the CB threshold V^^\ the interme- 
diate voltage region, and the large voltage region. Several static trajectories exist near the 
threshold, and a crossover from a static to a dynamic trajectory occurs when the bias voltage 
is set in the intermediate region. As the bias voltage increases, trajectories are again static 
and linear. The I-V property thus approaches to Ohmic behavior in the large voltage region. 

The I-V property near the CB threshold is mainly characterized by the value of the CB 
threshold. The threshold is determined by trajectory of electrons, and is thus sensitive to 
conditions such as the array size and surface disorder. The I-V property is approximately 
described as 

i-(v- (1) 

where V is the bias voltage. In the intermediate region, the I-V property also exhibits non- 
linear behavior described as eq. (1). The value of ( for several systems has been determined 
from both experiments and simulations. For example, an experimental study shows that an 
array of normal metal islands has ( — 1.36 ± 0.1 for a one-dimensional (ID) array and 
C = 1.8 ± 0.16 for a two-dimensional (2D) square array;^^ in addition, other experimental 
studies show that a metal nanocrystal has ( — 2.25 ±0.1 for a 2D triangle array,^''^ and that 
a gold nanocrystal has ( — 2.7 to 3.0 for a 3D array.^^^ Experiments of colloidal deposition 
show C = 2.1 for 2D and ( — 3.5 for 3D.^°^ For numerical simulations, MC calculations 
show that C = 10 for linear arrays and ( = 2.0 for square arrays,^^^ and MD calculations 
show that = 1.94 ± 0.15 for square array s.^^^ In a theoretical study, a mean-field analysis 
suggested that C = 2 for a 2D array.^^^ By analyzing surface evolution on arrays with charge 
disorder, = 5/3 is analytically predicted. ^'^ Here, we should emphasize that has been dis- 
cussed in relation to the array configuration and array dimension so far; the size dependence 
has not been taken into account. Meanwhile, some previous studies^' ^^'^"^^ mention that the 
exponent ( monotonically increases with increasing the lateral size from a ID linear array to 
a 2D square array. However, they have qualitatively focused on only the arrays in which the 
lateral size is less than the longitudinal one and expected the exponent to be constant for large 
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Fig. 1. The present configurations: (a) a simple lattice (SL), (b) a line-type triangular lattice (TL-1), and (c) a 
zigzag-type triangular lattice (TL-z). Each circle indicates a Coulomb island. This SL contains N = 4:8 islands 
{N^ = 8, Ny = 6, and S = 0.75), this TL-1 does N = 80 islands {N^ - 8, Ny = 10, and 5 = 1.25), and this 
TL-z does N = 72 islands (Nx = 12, A^j^ = 6 and S = 0.5). The islands are sandwiched between the positive 
(left) and negative (right) electrodes, and each solid line between island-island or island-electrode represents the 
tunneling junction. Each island touches the gate electrode (not displayed in the figure). 



lateral size. 

Middleton and Wingreen (MW) explicitly introduced offset charge distribution in their 
model. The charge disorder originates from the surface impurity. In their model, Bascones et 
al. have discussed the asymptotic I-V property of ID arrays in the large voltage region.^^^ It 
converges to the Ohm's law at the large voltage limit. In addition, they showed the presence 
of the offset voltage Kffset, and analytically expressed it in short-limit of the interaction range. 

In this paper, we carry out MC simulations to study the size dependence of the I-V 
property for configurations such as a simple lattice and a triangular lattice. We employ the 
model proposed by MW.^'^ Based on their model, we extend it to size dependence. Our main 
results are the following: (i) the average CB threshold , (ii) the power-law exponent ( 
in the intermediate voltage region, and (iii) an asymptotic I-V curve of 2D arrays with the 
first-order perturbation of e. 

This paper is organized as follows. In § 2, we briefly describe the present configurations, 
simulation model, and numerical conditions. In § 3, we first express the size dependence 
of the average CB threshold for simple configurations, and then, the size dependences of 
the exponent ( is shown for several configurations. In addition, we express the asymptotic 
I-V property obtained analytically for simple configurations, and then we compare it with 
simulation and experimental results to verify the asymptotic relation. In § 4, we summarize 
our results. 
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2. Method 

2.1 Structure 

The simplest single-electron transistor consists of two tunnel junctions that connect the 
source and drain electrode, respectively, and the region sandwiched between the tunnel junc- 
tions also connects to the gate electrode through the gate capacitor. The sandwiched region, 
known as the Coulomb island, can be regarded as a place where charge accumulates. We 
consider arrays constructed of Coulomb islands between positive and negative electrodes. A 
series of tunneling processes cause electrons to flow in the arrays. Each island also connects 
to the gate electrode with the gate capacitor. The i-th island has charge Qi and potential $i. 
The charge Qi contains both an integer multiple of the elementally charge ne (where n de- 
notes an integer and e the elementary charge) and offset charge — e/2 < qi < e/2 due to the 
impurity.^^' In simulations, the offset charges are set by uniform random numbers and remain 
constant over time. 

Three configurations are considered (Fig. 1). The first configuration is a simple lattice 
(SL), and the remaining two configurations involve different directions of a triangular lattice: 
a line-type triangular lattice (TL-1) and a zigzag-type triangular lattice (TL-z). We set x- 
and y-directions as shown in Fig. 1, namely x-direction corresponds to longitudinal direction 
and y-direction does to lateral one. The SL configuration is characterized by the number 
of horizontal islands and vertical islands Ny. Thus, the total number of islands is N — 
N^xNy. Both the TL-1 and TL-z configurations are also characterized by and A^^^, and the 
total number of islands is N = N^x Ny. Although we use A'^ and Ny in every configurations, 
the method of setting them is underspecified for triangular lattices. Here, we define A^^ and 
Ny as described in the caption of Fig. 1, and the aspect ratio S is defined as the lateral size 
over the longitudinal size; i.e., 6 = Ny/N^. 

To calculate the total energy of the array, it is useful to consider a configuration matrix. 
The physical configuration of the lattices uniquely determines the configuration matrix whose 
element is represented as 



where Cy denotes the tunneling capacitance between the i-th and j-th islands, the tunnel- 
ing capacitance between the i-th island and the electrode n G {+,—,g}, and 5ij the Kronecker 
delta. Note that the symbols — , and g indicate the positive, negative, and gate electrodes, 
respectively. If the i-th island does not connect to the j-th one, then Cy is set to 0. Electrons 




(2) 
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move through the network of islands, while the islands themselves do not move. Therefore, 
the configuration matrix Mjj remains constant over time in this study. 



2.2 Model 

We briefly summarize the time-evolution procedure used in the MC method^*^*^^^ in this 
subsection. 

The system evolves to decrease the total electrostatic energy E, whose derivation is sum- 
marized in Appendix A. In MC simulations, the electrons are virtually moved for each possi- 
ble tunneling event. We can calculate the energy change A£^n/^m' at n' m' (see Appendix 
B), where {n', m'} e {1, 2, ... , A^, +, -}. The tunneling rate T^i 

^m' n' — )■ m' is 

calculated as^^^ 

p ^ ^ ^ 1 -Ai?n'-^m' 

" e^i^t^n/^iji, 1 - cxp [/\En,^^,/kBT] ' 

where i?t,n'^m' denotes the tunneling resistance at n' — > m'. It is assumed that the tunnel 
resistance between an island and the gate electrode is infinity in these calculations, i.e., the 
electrons cannot move between an island and the gate electrode. The tunneling rate is derived 
by assuming that the tunneling events occur independently. The resistances i?t,n'-*^m' depend 
on the configuration of the array in general, but we regard them to be a constant i?t, where 

Assuming Poisson distribution, the probability distribution that a tunneling from n' to m' 
occurs at time lag t is represented as 

/n'^m' {t) = r„'^m' exp [-V^,_^^d] . (4) 

Then, the cumulative distribution is equivalent to the distribution that a tunneling from n' to 
m' occurs during time lag t — to, represented as 



F„/^.in'(i) = 1 - exp 



Jto 



(5) 



where to denotes the time when the last tunneling event from n' to m' occurred. The energy 
changes depend on time, and hence the tunneling rate also does. A uniform random number 
X = x{t) in [0, 1] is introduced and the cumulative distribution is set as F^i^^r{t) = x. Note 
that the random number x is updated only when a tunneling event from n' to m' occurs. We 
thus obtain the time interval between the last tunneling event from n' to m' and the next one 
as 

A. - l0g(l -X)- Ef="o^ rn/-^m'(^k)Atk 

iXtn'^m' = :^ J-. , iP) 

t n'^m'^lj 
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where K denotes the number of tunneUng events within the entire array during t — to. Note 
that rn/_^iii' remains constant over time At^, and Atn'-^m' is not equal tot — to in general. We 
use the smallest Atn'^m' for the time evolution increments. 

2.3 Simulation 

Below the CB threshold VJj^ , the tunneling interval At is infinity for each path. There- 
fore, we can determine the threshold voltage above which At is finite in the steady state. The 
current along the path n' m' can be calculated as 

In simulations, the current can be calculated along all path, but we focus only on those paths 
that neighbor the positive or negative electrodes, represented as 




i 



where the sigma with the prime denotes summation over only those islands that neighbor 
the positive or negative electrodes. In the steady state, /posi and /nega are the same because of 
Kirchhoff 's current law. Hence, we demonstrate only /posi as the current / in the remaining 
sections. 

The voltages of the negative and gate electrodes are fixed at = $g = 0, and the voltage 
of the positive electrode is thus adjusted as the control parameter, i.e., the bias voltage V 
is equivalent to the voltage of the positive electrode $+. The initial condition was Qi = 
and $+ = 0. The voltage was incremented by A$+, and before we sampled the physical 
variables at each voltage, we waited sufficiently long for the system to return to the steady 
state. Note that this waiting time depends on system conditions, such as the configuration and 
A$+. 

We assume that the system is at zero temperature. The capacitance is set at C = 10~^Cg, 
Note that the ratio e :— C/Cg corresponds to the interaction range.^^^ The increment voltage 
A$+ is 10~^; therefore the threshold voltage has an uncertainty of the order 10~^. Finally, 
the charge is scaled by e, the capacitance by Cg, the time by RiCg, the current by e/ RiCg, the 
potential by e/C^, and the energy by /Cg. 

3. Results and Discussion 

3.1 Size dependences of the average threshold 

We first investigate the size dependences of the average CB threshold for SL. Figure 2 
shows the average CB threshold V'th*^^^ as a function of the aspect ratio 5 for several longitudi- 
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Fig. 2. (Color online) Log-log plot of average CB threshold V^l^^ versus the aspect ratio S for several longi- 
tudinal sizes; ~ 1 (circle), 2 (square), 5 (triangle), and 10 (diamond). The red solid line represents eq. (12), 
and the red dashed line represents power-law fitting. The data in 5 > 5 are used to obtain the fitting parameters. 





Fig. 3. (Color online) The coefficient c of eq. (13) versus a length of a side of square arrays. A filled circle 
represents a simulation result. The red dashed line is a line with slope 0.338. 




Fig. 4. (Color online) The exponent 7 of eq. (13) as a function of (i.e., S at fixed Ny). A filled square 
represents a simulation result. The dashed line represents eq. (14). 
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nal size N^. Each point on the figure is derived from the average of at least 50 different initial 
distributions of the offset charges. Our results of ID arrays (i.e., Ny = 1) are in agreement 
with the previous results by MW.^^^ 

The average threshold is analytically represented as 

vr\N., N,) = r ^ ■ • ■ r ^v-r'({«,}). (9) 



J-e/2 6 J -e/2 ^ 

For the simplest case = Ny = 1 and £ <^ 1, the CB threshold voltage as a function of the 
initial charge q is obtained as 

namely, V^*^^^ of a single Coulomb island is proportional to the initial offset charge. For 
Nx — 1 (i.e., N — Ny) and £ <^ 1, the threshold is dominated by the smallest initial offset 
charge. The average threshold reduces to 

V,^^\N^ = l,Ny) 



J-1/2 7-1/2 7-1/2 \ ^/ 

(11) 

where {q{} denote the reordered dimensionless charges: 1/2 > q[ > q'2 >■■■> q'j^ > —1/2. 
The average threshold for AT^ = 1 is thus obtained as 

Vr(N. = l,N,) = ^^^^, (12) 

and this well describes the simulation result as shown in Fig. 2. For Nx > 1, it is difficult to 
derive the average threshold because electron meandering plays an important role just above 
VJjj . Nevertheless, we find that the average CB threshold for large 5 can be described by a 
power law 

vi''^\Nx,Ny)^c5-\ (13) 

In fact, eq. (12) for large Ny implies the power-law decay as eq. (13). 

The proportionality coefficient c of eq. (13) indicates the value of the average threshold 
for square arrays (i.e., Nx — Ny). Figure 3 shows the coefficient c as a function of Ngq that is 
a length of a side of square arrays. MW have reported V^^^^ ~ 0.338A^sq for square arrays,^^^ 
and the line is plotted in Fig. 3. The simulation results deviate from the line in small Ngq 
region. This is because we cannot regard the array with A^j; = A^j, = 1 as a square array, 
namely, N^ is too small to regard arrays as square in that region. In addition, the simulation 
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Fig. 5. (Color online) Average current is shown for SL with Nx = 40. (a) I-V plot for the aspect ratio 
6 = 0.025, 0.125, 0.25, 0.5, 1, and 2. (b) Current as a function ofV - V^^^^ for the same values of 6. 



results do not satisfy the power law near 5 ~ 1 at small A^^;. 

The power-law exponent 7 of eq. (13) is shown in Fig. 4. The relation 7 ~ 1 for A^^. = 1 
is evident from eq. (12). The exponent 7 is thus interpreted as a sensitivity of the threshold 
depending on 5 comparing to arrays with A'^,. = 1. As increases, the increment of the 
aspect ratio decreases even for the same increment of Ny. Therefore, it can be expected that 
7 is a monotonically decreasing function of N^- In fact, as shown in Fig. 4, the exponent 7 is 
inversely proportional to N^, namely, 

^{N,)=N^-\ (14) 

These results will be a hint to understand the size dependences of the CB threshold. Li addi- 
tion, it is interesting to observe the power law decay in experiments. 

3.2 Logarithmic increase of the power-law exponent ( 

We next investigate the exponent ^ in eq. (1) as a function of the aspect ratio 5. Figure 
5 shows the averaged I-V property for several 5, and each curve results from the average 
of at least 30 data sets. The nonlinear behavior is more visible in Fig. 5 (b). Further, Fig. 6 
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Fig. 6. (Color online) Semi-log plot of the exponent ( for the simple lattice as a function of the aspect ratio 
S with Nx = 40. This plot is extracted from the average I-V relation shown in Fig. 5 (b). The red dashed line 
shows the fitting fine determined by eqs. (15) and (16) with C*"*^ = 2.08 and C*""'^ = 1.05. 

shows the exponent C as a function of 6 with A^^ = 40. The exponents are obtained by fitting 
to the average I-V property in the intermediate region defined as 10°'^ < V — V^j^^^^ < 10. 
The fitting range is selected to prevent artificiality from being included into the value of the 
exponents. The aspect ratio dependence of appears to be approximately represented by 



with fitting parameters C^^''^ and b = b(Nx). The C^^''^ parameter denotes the exponent of the 
square array, and (^^^^ ~ 2.08 approximately agrees with the previous result. Using the b 
parameter, the exponent ^('i"^) for a ID simple array, i.e., Ny = 1, is represented as 



and ^('1"*=) ~ 1.05 is obtained as shown in Fig. 6. MW have reported that, for arbitrary 
N^, the exponents of linear and square arrays are C — 1 ^^id ^ ~ 2, respectively.^^' The 
same logarithmic increase is thus expected to appear in different longitudinal size (i.e., A^^.)^ 
while we only show single longitudinal size (i.e., A^^; = 40) in Fig. 6. Moreover, although 
our simulations are only performed for finite 6 in each configuration, it is expected that the 
exponent diverges logarithmically at such large values of 5. 

We also focus on the triangular lattices. Figures 7 (a) and (b) show the exponent as a 
function of the aspect ratio S for TL-1 and TL-z, respectively. Note that, similar to SL, the 
I-V properties are averaged over at least 30 data sets and the exponent is extracted from the 
averaged results. Here, the intermediate region is defined as 10°-^ < V — VJh™^ < 10°'^, 



C = C^^^) + bhg,,S 



(15) 



^(line)^^(sq)_^l^g^^^^ 



(16) 
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Fig. 7. (Color online) Semi-log plot of the exponent ( for (a) TL-1 and (b) TL-z as a function of the aspect ratio 
6 with Nx = 20. The red dashed lines show the fitting lines determined by eqs. (15) and (16) with (^'^^''^ — 1.86 
and C^""") = 0.884 for TL-1 and with C^"*^ ^ 1.96 and C^'^^ = 1.06 for TL-z. 



and the fitting is done in that range. Similar to the exponent of SL, the exponents ( of TL-1 
and TL-z increase logarithmically as determined by eqs. (15) and (16). Parthasarathy et al.'"^ 
experimentally showed ( ~ 2.25 for a well-ordered triangular array of gold nanocrystals 
with an array size of A^^^^ = 30 to 90 and Ny ~ 270. In the range 3 < 5 < 9, our simulation 
results of both TL-1 and TL-z show ( ~ 2.25. Our results propose that we should pay attention 
to the aspect ratio as well as the array configuration and array dimension when discussing the 
exponent (. 

We should recall that the universality of the logarithmic increase requires careful atten- 
tion. The above results are obtained in locally coupled CB, i.e., small e. The behavior of C for 
large e may differ from that for small e, because the interaction among electrons ranges over 
the entire array in large e systems. The e dependences are still an open question. 
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3.3 Analytical asymptotic equation at large bias voltages 

We finally discuss the I-V property for large values of the bias voltage y = — 
As mentioned above, Bascones et al. have derived the asymptotic I-V property of ID arrays 
with the offset voltage Kffset Their offset voltage Kffset is expressed in the limit of £ — )■ 0.^^^ 
We extend their study to 2D simple configurations. In addition, our extended expression of 
Kffset contains the first-order perturbation of e. 

Figure 8 shows the I-V plot for several sizes of the SL configuration. All results exhibit 
linear behavior at large bias voltage limit, and the asymptotic curve can be obtained as fol- 
lows. The energy changes of the ID array that contains islands reduce to (see Appendix 
B) 

2 

AE^^+ = e + ^/'''•^) - e$+ + yMfiS (17a) 

A£;„+i_„ = e(K+i + Kl?-K-K^'''')) 

+ ^ (M-„^ + M„-\,„+i - 2M-„Vi) , (17b) 

2 

with n = 1, 2, . . . , A^^ — 1. At large bias voltages, the current from the negative to the positive 
electrode is neglected. Thus, the all energy changes should be always the same in the large 
voltage region. We can obtain AE' := AEi^+ — AE2^i — ■ ■ ■ — AE ^jv^, as 

AE^--^{V-Vomei). (18) 
Here, the above representation contains an offset voltage Kffset defined as 

The offset voltage Kffset differs from the CB threshold V^ . For e <^ I, Kffset can be ana- 
lytically derived within the first-order perturbation of e (see Appendix C) as 

K.set - ^ ^1 j . (20) 

The energy changes cannot be the same below the offset voltage Kffset- In contrast, above 
the offset voltage Offset, electrons can move from the negative to positive electrode for any 
offset charge distributions, i.e., Kffset is the maximum of V^^^ when e <S 1.^^^ However, near 
T4ffset> the influence of the offset charge distribution cannot be neglected because the charge of 
electrons is discrete. With increasing the bias voltage, the energy changes become sufficiently 
large to neglect this discreteness. Therefore, all energy changes in the large voltage region can 
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be always regarded as the same. 

Using eq. (7), the asymptotic current is obtained as 



j(asy, ID) 



V -Vr 



offset 



1 — exp 



e{V-V, 



offset J 



(iV,. + l)kBT 



n -1 



(21) 



Note that it is independent of the gate voltage. The simple array can be easily extended to 
higher-dimensional arrays. For the 2D array in which the number of islands is A^^, x Ny, the 
lateral current (i.e., the ^/-direction current in Fig. 1) can be neglected at large bias voltages. 
This 2D array can be assumed to be composed of the isolated A^j, ID arrays that consist of 
Nrr islands. Hence, we obtain 



j(asy, 2D) 

Ny iV - Kifset) 

(iV. + l)i?t 



exp 



e{V-V, 



offset J 



(iV^ + l)kBT 



(22) 



The above equation does not hold when ksT — )■ oo, because the assumption that the current 
from the negative to positive electrode is negligible is no longer correct at these high tem- 
peratures. In contrast, in finite temperature, eq. (22) represents that the asymptotic equation 
converges in the limit of V^/Kffset — oo to the Ohmic behavior V ^ R^I with the combined 
tunneling-resistance 



Rr 



(23) 



Namely, the combined tunneling-resistance is inversely proportional to the aspect ratio 5; 
RjRt ~ S-^ at large N^. 
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Fig. 9. (Color online) Current- voltage plot for several sizes {Nx,Ny) of the SL configuration at T = 
for the results shown in Fig. 8; the asymptotic line calculated from eq. (22) is shown by the black dashed 
line. The horizontal axis denotes V — Kffset and the vertical axis denotes R^I with the combined resistance 
i?c = {Nx + l)Ri/Ny. All results for large V collapse to the asymptotic line. 



As shown in Fig. 9, which plots R^I as function of ^ — Kffset, the asymptotic line calcu- 
lated from eq. (22) completely describes simulation results of arbitrary S at large bias voltages 
(roughly, V^/Kffset > 2). Near the offset voltage Kffset, the simulation results deviate from the 
asymptotic line because each energy change is different from the others. As mentioned above, 
this originates from discreteness of charge. In fact, as A^^^ decreases (i.e., the number of the 
energy changes decreases), the results collapse to the asymptotic line at smaller V — Kffset- 

We next compare the asymptotic result with the experimental result^^ measured by Kurdak 
et al. Because the physical parameters are known, this experimental result is suitable for 
testing the asymptotic equation. Figure 10 shows the experimental (sample A in reference^^) 
and the asymptotic results for the following experimental conditions from reference: = 
Ny = 40 (i.e., 6 = 1), Cg = 1.38 fF, C = 0.25 fF, = 810 kQ, and T = 20 mK. 
The temperature is sufficiently small for neglecting the exponential dependence in eq. (22), 
and we neglect the square of e = 0.0181. The asymptotic results typically describes the 
experimental result, as shown in Fig. 10. The asymptotic fitting parameters lead to -Rt = 789 
\dl and C = 0.0128 fF. The resistance -Rt is in good agreement with the value in reference.^^ 
On the other hand, the capacitance C evaluated from the asymptotic equation is smaller than 
that estimated by Kurdak et al. Several reasons can be considered. First, our expression of 
Kffset neglects higher terms of e and one might consider the higher-term effects. Second, the 
capacitance ratio e is sensitive to Kftset- As shown in Fig. 10, the order of C obtained from 
fitting is different from that estimated by Kurdak et al., while the difference of Kffset is 2.4 mV. 
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Fig. 10. (Color online) Current-voltage property for the experimental result extracted from figure 2 (sample 
A) in reference^' (black solid curve) and the analytical result calculated from the following formula and values 
given in reference^* (blue dotted line): / — a{V ~ Kffset) with a = 1.20 /iA/V and Kffset = 2.17 mV. In 
addition, the asymptotic fit for large voltages (a = 1.24 fiAfV and Kffset — 4.57 mV) is shown (red dashed 
line). Note that the experimental result is digitized from reference.^' 



In addition, it is difficult to experimentally evaluate the value of C in general. In fact, Kurdak 
et al. stated "our knowledge of C is less precise," and they estimated the value of C from 
the specific capacitance. Instead, the asymptotic equation may allow us to approximately 
evaluate the configuration of the array and some physical variables (N^, Ny, Ri, C and Cg) in 
experiments from observations of the offset voltage Kffset and the asymptotical slope at large 
voltages. 

4. Summary 

We conducted MC simulations to investigate the I-V properties of CB arrays. To under- 
stand the I-V property, our strategy was dividing it into three regions characterized by the 
path flow of electrons, and we pay attention to the size (i.e., aspect ratio 5) dependence. Our 
main results were (i) power-law behavior of the average CB threshold V^^\ (ii) the power- 
law exponent C, in the intermediate voltage region, and (iii) an asymptotic I-V curve at large 
voltages. 

We derived an analytical relationship [eq. (12)] for the average CB threshold for = 1, 
and we found that the average CB threshold obeys a power-law decay as a function of the 
aspect ratio S. The coefficient c is in agreement with the previous study by MW.^^^ In addition, 
its power-law exponent 7 is inversely proportional to the longitudinal size (i.e., 5 at fixed 
Nx) [eq. (14)]. It is difficult to obtain the analytical form for > 1 because trajectories 
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of electrons meander. Nevertheless, our analytical and simulation results provide a hint for 
further development of the study. 

The size dependences of the exponent ( were shown for different array configurations 
such as SL, TL-1, and TL-z. The exponent in arrays of large 5 was considered to be constant 
so far. However, we revealed that ( logarithmically increases as 6 increases for both the simple 
and triangular lattices. Namely, in addition to the array configuration and array dimension, 
the aspect ratio 5 is a significant variable for discussing the exponent (. 

We extended the asymptotic equation for ID arrays without interaction^^^ to 2D arrays 
with first-order perturbation of the interaction range e [eq. (22) with eq. (20)]. At suffi- 
cient large voltages, the equation adequately describes the Ohmic behavior and the combined 
tunneling-resistance Rc is inversely proportional to 5. The offset voltage Kffset^ included in 
the asymptotic equation, differs from the CB threshold V^^^\ Instead, the offset voltage can 
be regarded as the maximum of the V^th*^^^ in the limit of e — 7> 0.^-'^ These asymptotic property 
well agrees with simulation and experimental results. Our extended equation allows to esti- 
mate physical values and array configuration which are experimentally hard to obtain. The 
asymptotic equations for other configurations are also expected to show similar results and to 
converges to the Ohm's law. The details will be discussed elsewhere. 
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Appendix A: Total Energy 

The charge of the i-th island is represented by 

gi = Yl (A-i) 

j M=+.-,9 

Assuming that the capacitances are nonzero only between neighboring island-island and 
island-electrode pairs, eq. (A l) reduces to 

j l^=+,-,9 

where Mjj denotes the matrix of capacitances defined by eq. (2). The capacitance Ca should 
be zero by definition. Equation (2) thus indicates that the diagonal elements Ma are the sum 
of all capacitances associated with an island, and the off-diagonal elements (i^^j) are the 
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negative of the capacitance between the i-th and the j-th islands. The potential $i is formally 
solved to obtain 

*i = I]Mr'^j+K^"^ (A-3) 

j 

where l^*^'^"'^ denotes the potential corresponding to the electrodes defined by 

j /*=+,->fl 

The total electrostatic energy of the system is equivalent to the sum of the work for storing 
charge Qi under potential $i in each island and the energy of the electrodes, represented as, 

where the last term denotes the energy of the electrodes and the charge at the electrodes 
/i G {+, — , (7}. Note that the interparticle electrostatic energy must not be double-counted. 

Appendix B: Energy Change 

Let us consider the tunneling of an electron whose charge is — e from the n-th to m-th 
island. The energy change is 

AE,^^ AE^Xra + ^Ei^-^l (B- 1) 

where AE^Xm and AE^^]a denote the energy changes with respect to the first and second 
terms of eq. (A-5), respectively. The charge changes to Q[ = Qi + e6in — eSim with the 
tunneling; therefore, 

i.j 
i 

+ - [M-,' + M-l-2M-2] 
= e[K-Kn] + ^[M^' + M^-2M^^], (B-2) 
where the trivial relationship M-^^ = M^^ is used and the effective potential is introduced as 

j j 



= e 

o2 
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Similarly, 

i i 

= e [K^*^*) - y^^"*)] . (B-4) 

Next, let us consider the tunneling from the n-th island to an electrode /j,. Similar to the 
above discussion, the energy change is represented as 

AE,^, := AE^X^ + AE^^l + AE^^^'^^^'l (6-5) 

Given that the charge change is Q[ = Qi + edin, the energy change is obtained as 

AE^X^ = lj2{Q, + e5,,)M,r\Q^ + e5i,,) 

= eK + yM„-i (B-6) 

and 

AEi':^l = eV}''''\ AS^!!,^*°'^'=) = -e$^. (B-7) 
In addition, the following equations hold: 

AE^^X^ = -eK.+ yM^, (B-8) 

AEj^^l = -eV^'^'\ (B-9) 

Appendix C: Offset Voltage Votiset 

The configuration matrix M^^^^ for a ID simple array is represented as 

r {l + 2e)C, i = j 

M^^'^^^i -£C, = l (C-1) 
[ otherwise 

with arbitrary e = C/Cg. The inverse elements M^-^^ ^ (i = 1, 2, • • • , A^^) and M^-^l ^ 
(j = 1, 2, • • • , A'a; — 1) are derived as 

Mjjj^i"^ = eAj_iA^_j_i/A^,, (C-3) 
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where A„ denotes the determinant of the configuration matrix for the ID simple array that 
contains n islands and Aq = 1. The offset voltage reduces to 

Kffset =^[Y1 ^i-l^AT.-i -^Y. ^i-l^iV.-i-l I . (C-4) 
\ i=l i=l / 

Note that the above representation holds for arbitrary s. 

We can approximately obtain the determinant A„ for £ ^ 1 as 

A„ = [{l + 2e + 0{e'))C,Y 

= (1 + 2ne + 0(£'))C/. (C-5) 
Substituting eq. (C-5) into eq. (C-4) leads to eq. (20). 
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